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1.0  SUMMARY 


The  global  advancement  and  proliferation  of  high  perfonnance  radio  frequency  (RF) 
commercial-off-the-shelf  (COTS)  technologies  has  created  significant  challenges  to  electronics 
intelligence  (ELINT)  and  Communication  Intelligence  (COMINT)  operations.  As  a 
consequence,  entirely  new  methods  of  signal  exploitation  must  be  developed  and  fielded  to  meet 
these  challenges.  This  work  focused  on  multi-sensor  algorithmic  approaches  for  the  geo-location 
of  LPI  signals  in  environments  with  significant  multipath  and  co-channel  signals.  While  many 
classical  techniques  exist,  new  methods  are  needed  that  will  allow  accurate  location  of  LPI 
emitters  despite  low  signal-to-noise  ratio  (SNR)  and  multipath  propagation.  Geo-location 
systems  based  on  time-difference-of-arrival  (TDOA)  and  frequency-difference-of-arrival 
(FDOA)  estimate  the  TDOA/FDOA  between  a  pair  of  signals  by  computing  their  cross¬ 
ambiguity  function  (CAF)  and  finding  the  location  of  the  peak  on  the  surface  of  this  2-D 
function.  However,  when  multipath  signals  are  received  at  each  platform  they  give  rise  to 
spurious  peaks  on  the  CAF  that  can  perturb  the  location  of  the  true  peak,  especially  if  the 
spurious  peaks  are  located  close  to  the  true  peak. 

This  classical  approach  uses  two  stages  to  estimate  the  signal  position.  In  the  first  stage, 
TDOA/FDOA  are  estimated  by  several  pairs  of  sensors  and  then  used  in  the  second  stage  to 
locate  the  emitter.  However,  this  two-stage  method  is  known  to  perform  poorly  in  low  SNR 
cases  and  especially  in  the  presence  of  multipath  and  multiple  co-channel  emitters.  To  address 
the  drawbacks  of  the  two-stage  method,  recent  work  has  proposed  a  single-stage  method  based 
on  TDOA/FDOA;  that  method  is  called  Direct  Position  Determination  (DPD).  However,  that 
improvement  comes  at  a  cost  of  significantly  more  computational  complexity,  and  unfortunately 
that  complexity  is  all  concentrated  at  one  computing  node,  unlike  in  the  classical  method  where 
the  computations  are  distributed  evenly  among  the  sensors.  As  originally  proposed,  the  single- 
stage  method  transmits  all  data  to  a  single  sensor,  which  then  forms  a  myriad  of  matrices  (one  for 
each  possible  emitter  location  grid  point)  and  computes  the  maximum  eigenvalue  of  each  one. 
The  largest  of  these  maximum  eigenvalue  then  indicates  the  emitter  location. 

We  developed  several  methods  to  ease  the  application  of  the  new  single-stage  method.  In 
particular,  we  show  how  to  reduce  the  load  of  data  computation  and  data  transmission  using 
distributed  data  computation  and  processing,  applying  data  compression  methods,  exploiting  the 
CAF  properties,  taking  advantage  of  CAF  relationships  in  the  sensor  network  and  exploiting 
some  kind  of  beneficial  approximations.  We  proposed  three  alternative  processing  schemes. 
The  Approximated  DPD  method  exploits  the  simplicity  of  Gershgorin’s  theorem  to 
approximately  compute  the  maximum  eigenvalue  without  the  high  cost  of  exactly  computing  it. 
This  enables  each  sensor  to  locally  make  its  best  estimate  of  the  location  based  on  that  data  it 
has.  These  locally-generated  estimates  are  then  transmitted  to  a  central  location  where  a  final 
decision  is  made.  In  Decentralized  DPD  method,  we  applied  the  distributed  computation  idea  to 
divide  the  mathematical  calculation  load  among  all  sensors  of  the  network.  In  this  method,  we 
don’t  use  any  other  approximation  more  than  data  compression  and  the  quality  of  this  method  is 
the  same  as  the  Original  DPD  for  reasonable  bit  rates.  Finally,  in  Semi-Optimal  Decentralized 
DPD,  we  used  the  same  idea  of  Decentralized  DPD  in  addition  to  exploiting  the  relationship 
between  different  CAFs  to  eliminate  the  data  redundancy  and  this  idea  leads  to  a  large  data 
transmission  reduction.  All  of  the  proposed  methods  allows  DPD  to  be  implemented  in  a 
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decentralized  manner  where  no  single  sensor  is  required  to  do  an  unfair  share  of  the 
computations,  yet  the  perfonnance  improvement  of  DPD  is  not  sacrificed. 

Although  the  DPD  method  provides  improved  location  performance  at  low  SNR  it  is  not  as 
robust  as  desired  in  the  presence  of  multipath  and/or  multiple  co-channel  emitters.  To  provide 
better  performance  in  these  challenging  real-world  scenarios  we  developed  a  one-stage 
TDOA/FDOA  localization  method  based  on  spatial  sparsity  of  emitters.  In  this  method,  we 
imagine  assigning  a  non-zero  number  to  each  one  of  the  grid  points  containing  an  emitter  and 
zero  to  the  rest  of  the  grid  points.  Thus,  the  vector  formed  from  these  numbers  will  be  a  sparse 
unknown  vector  that  we  aim  to  estimate.  Since  each  element  of  this  vector  corresponds  to  one 
grid  point  in  the  (x,y)  plane,  we  can  estimate  the  location  of  emitters  by  extracting  the  position  of 
non-zero  elements  of  the  sparsest  vector  that  satisfies  the  TDOA/FDOA  relationship  between 
transmitted  signals  and  received  signals.  The  proposed  sparsity-based  method  has  better 
performance  (especially  in  multi-path  and  multi-emitter  cases)  compared  to  direct  position 
detennination  (DPD)  and  two-stage  Classic  localization  methods. 


2.0  INTRODUCTION 

The  global  advancement  and  proliferation  of  high  performance  radio  frequency  (RF)  com¬ 
mercial-off-the-shelf  (COTS)  technologies  has  created  significant  challenges  to  electronics 
intelligence  (ELINT)  operations.  Typically  ELINT  refers  to  RADAR  signals,  however  the  same 
issues  have  also  faced  Communication  Intelligence  (COMINT).  For  example,  technology 
advances  such  as  digital  arbitrary  wavefonn  generators  (DAWGs),  solid  state  transmitters,  active 
electronically  scanned  arrays  (AESAs),  high  speed  analog-to-digital  converters  (ADCs),  and 
spread  spectrum  low-probability-of-intercept  (LPI)  waveforms,  have  given  rise  to  a  new  class  of 
signals  that  are  extremely  difficult  to  detect,  de-interleave  and  characterize.  While  ELINT 
receivers  have  increased  their  operating  bandwidths,  that  alone  is  insufficient  to  guarantee 
successful  intercept  as  this  is  a  necessary  but  not  sufficient  condition.  As  a  consequence,  entirely 
new  methods  of  signal  exploitation  must  be  developed  and  fielded  to  meet  these  challenges. 

This  work  focused  on  multi-sensor  algorithmic  approaches  for  the  geo-location  of  LPI 
signals  in  environments  with  significant  multipath  and  co-channel  signals.  While  many  classical 
techniques  exist,  new  methods  are  needed  that  will  allow  accurate  location  of  LPI  emitters 
despite  low  signal-to-noise  ratio  (SNR)  and  multipath  propagation.  Geo-location  systems  based 
on  TDOA/FDOA  estimate  the  TDOA/FDOA  between  a  pair  of  signals  by  computing  their  cross¬ 
ambiguity  function  (CAF)  and  finding  the  location  of  the  peak  on  the  surface  of  this  2-D 
function.  However,  when  multipath  signals  are  received  at  each  platform  they  give  rise  to 
spurious  peaks  on  the  CAF  that  can  perturb  the  location  of  the  true  peak,  especially  if  the 
spurious  peaks  are  located  close  to  the  true  peak. 


3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

One  of  the  most  accurate  and  common  methods  for  passive  radio  signal  geolocation  is  based 
on  TDOA/FDOA  estimation.  The  classical  approach  to  this  method  uses  two  stages  to  estimate 
the  signal  position.  In  the  first  stage,  frequency-difference-of-arrival  (FDOA)  and  time- 
difference-of-arrival  (TDOA)  are  estimated  from  the  cross-correlation  of  signals  received  by 
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several  pairs  of  sensors  [1];  this  is  done  by  computing  the  cross  ambiguity  function  (CAF)  [2] 
and  finding  the  peak  of  its  magnitude  surface  [1],  [2].  In  [4]  and  [5],  a  Fisher  Infonnation  based 
data  compression  method  has  been  suggested  to  reduce  the  amount  of  data  transmission  and 
improve  the  communication  perfonnance  between  each  pair  of  sensors.  In  the  second  stage  of 
the  classic  method,  the  TDOA/FDOA  estimates  are  used  in  statistical  processing  to  locate  the 
emitter  [3]. 

Suppose  that  the  lowpass  equivalent  (LPE)  model  of  the  received  signal  is: 


sr(t)  =  aeJ°vs(t  -rd)  +  v(t) ,  (1) 

where  s(t )  is  the  LPE  of  the  transmitted  signal,  Wd  is  the  Doppler,  Zd  is  the  delay  for  the  received 
signal,  or  is  a  complex  number  and  v(t)  is  the  LPE  of  the  noise  [11].  Now,  suppose  that  two 
sensors  R1  and  R2  receive  the  LPE  signals  srl(t)  and  sr2(t) ,  respectively.  Stein  [1]  showed  that 

the  maximum  likelihood  (ML)  estimate  for  TDOA  and  FDOA  can  be  obtained  by  finding  the 
peak  of  the  magnitude  of  the  CAF: 


r  +0°  *  ■  * 

CAFn(z,co)  =  [  srX(t)sr2(t-z)e'  dt ,  (2) 

J  -oo 

which  measures  the  correlation  between  srl(t)  and  a  Doppler-shifted  by  o>  and  delayed  by  z 
version  of  sr2(t) .  The  accuracy  of  the  first  stage  is  governed  by  the  Cramer-Rao  lower  bounds 
(CRLB)  for  TDOA  and  FDOA.  Stein  [2],  Wax  [20],  Fowler  and  Hu  [21]  and  Yeredor  and  Angel 
[22]  derived  formulas  for  the  CRLB  on  TDOA  and  FDOA.  Yereder  [22]  have  derived  the  CRLB 
for  general  unknown  deterministic  signals  and  the  simulation  results  show  that  his  approach 
obtains  more  accurate  results  compared  to  other  CRLB  fonnulas. 

Recently,  some  new  methods  based  on  TDOA/FDOA  emitter  location  have  been  proposed 
that  estimate  the  emitter  location  in  one  stage  without  extracting  the  TDOA/FDOA  in  a  separate 
stage.  The  goal  of  these  methods  is  to  improve  the  overall  accuracy  of  the  emitter  location 
estimate  especially  in  low  SNR  cases.  Weiss  and  Amar  [7],  [8],  [9]  showed  that  the  two-stage 
methods  are  not  necessarily  optimal  because  in  the  first  stage  of  these  methods,  the  TDOA  and 
FDOA  estimates  are  obtained  by  ignoring  the  fact  that  all  measurements  should  be  consistent 
with  a  single  emitter  location.  In  other  words,  we  can  say  that  each  TDOA/FDOA  estimation  is 
optimal  in  the  first  stage.  Also  in  the  second  stage,  the  location  estimation  is  also  optimal  based 
on  TDOA/FDOA’ s  obtained  from  the  first  stage.  But,  it  does  not  necessarily  mean  that  the  whole 
two-stage  method  is  optimal. 

In  related  work,  Kay  and  Vankayalapati  [19]  developed  the  generalized  likelihood  ratio 
(GLR)  detector  based  on  the  received  signals  from  all  sensors  and  the  DPD  location  result 
appears  as  the  ML  estimate  used  in  the  GLR.  This  shows  another  advantage  of  DPD  over  the 
classical  two-stage  method:  the  classical  method  can’t  make  use  of  the  data  from  a  CAF  whose 
peak  is  undetectable  due  to  low  SNR  -  yet  the  DPD  method  can. 

However,  there  are  some  issues  that  need  to  be  addressed.  Namely,  the  implementation  of 
the  DPD  method  has  many  challenges,  such  as  how  to  distribute  the  computation  across  the 
participating  sensors  and  how  to  efficiently  communicate  the  necessary  data  between  the  sensors. 
Furthennore,  the  performance  of  the  DPD  method  in  the  presence  of  multipath  and  multiple  co¬ 
channel  emitters  was  previously  unknown.  This  report  addresses  the  distributed 
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computation/communication  issues  of  DPD  and  also  demonstrates  that  DPD  perfonnance  can 
suffer  in  the  presence  of  multipath  and  multiple  co-channel  emitters.  To  address  these 
shortcomings  we  have  extended  the  general  single-stage  idea  to  be  able  exploit  the  sparsity  of 
emitters  in  the  geometrical  region. 

4.0  RESULTS  AND  DISCUSSION 

4.1  Communication  and  Distributed  Computation  for  Single-Stage  TDOA/FDOA 
Location 

In  the  DPD  method,  we  need  all  the  received  signals  together  at  a  single  point  to  start  the 
location  estimation  processing.  Consequently,  all  sensors  have  to  transmit  their  received  signals 
to  a  common  site,  which  usually  is  one  of  the  sensors  so  we  will  refer  to  this  as  the  common 
sensor.  The  common  sensor  then  uses  the  received  signals  to  form  a  series  of  matrices  (one  for 
each  point  on  an  x-y  location  grid)  and  computes  the  maximum  eigenvalue  of  each  of  these 
matrices;  the  location  estimate  is  the  grid  point  that  produced  the  largest  of  these  maximum 
eigenvalues. 

As  mentioned  above,  the  one-stage  DPD  method  achieves  more  accurate  results  compared  to 
classic  two-stage  methods.  However,  in  the  published  papers,  the  authors  did  not  address  the 
issues  of  computation  and  data  transmission  for  DPD  and  there  are  some  difficulties  that  may 
limit  DPD  applications  in  practice.  The  first  problem  is  the  large  amount  of  computations  that 
are  to  be  done  by  only  the  common  sensor.  As  mentioned  above,  in  DPD  method  we  need  all 
received  signals  together  to  start  the  estimation  process.  Thus,  all  sensors  should  send  their 
received  signals  to  one  common  sensor  to  start  the  estimation  process.  The  common  sensor  will 
do  all  mathematical  computations  having  all  received  signals.  This  leads  to  a  large  computational 
load  on  only  one  point  in  the  network  and  no  computational  load  on  other  members  of  the  sensor 
network,  which  requires  any  one  sensor  in  the  system  to  be  computationally  capable  of  doing  the 
complete  set  of  computations  needed  to  locate  an  emitter.  In  scenarios  where  the  amount  of 
computational  capabilities  any  one  sensor  possesses  is  limited  this  centralized  approach  is  not 
desirable. 

The  second  problem  is  the  large  amount  data  transmission  to  one  single  point.  In  other 
words,  large  bandwidth  data  links  are  required  to  transfer  all  received  signals  to  the  common 
sensor  and  it  leads  to  have  a  bottle-neck  at  the  common  sensor  channel. 

The  third  problem  is  the  high  dependence  of  the  whole  network  on  the  common  sensor.  In 
this  scenario,  if  we  lose  the  common  sensor  during  computations,  we  will  lose  everything.  In 
other  words,  it  is  not  desirable  to  rely  on  only  one  point  in  the  sensor  network  for  computations 
and  data  collection,  because  if  we  lose  that  sensor  for  any  reason,  then  we  will  lose  all 
intennediate  and  final  results. 

In  this  paper,  we  develop  some  methods  to  increase  the  flexibility  and  feasibility  and 
improve  the  performance  of  one-stage  geolocation  methods.  We  use  distributed  data 
compression  to  reduce  the  amount  of  data  transmission  in  suggested  methods  and  also  we  use 
distributed  computations  in  the  sensor  network  to  reduce  the  mathematical  computational  load  on 
the  common  sensor.  Comparisons  will  be  made  to  detennine  the  advantages  and  disadvantages 
of  each  method  in  terms  of  estimation  accuracy,  computation  load,  transmission  load  and 
reliability. 
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4.1.1  Background 

In  this  section  we  provide  some  more  details  about  the  DPD  single-stage  localization  method  and 
its  formulation  [9].  Suppose  that  there  are  L  moving  sensors  in  a  sensor  network  receiving  the 
transmitted  signal  in  one  short  snapshot.  The  complex  signal  observed  by  the  /th  sensor  is 

Srl  (t)  =  ats(t  -Tl)ej2nf,t  +  w,(t ) 

where  s{t )  is  the  transmitted  signal,  a,  is  an  unknown  complex  path  attenuation,  ,  f  is  the 
Doppler  shift,  r,  is  the  signal  delay  and  w,  (t  )  is  a  white,  zero  mean,  complex  Gaussian.  Assume 
that  each  sensor  collects  A  time  samples  sampled  with  sampling  frequency  Fs  =  l/Ts  .Then,  we 
have 


srl  =  a ,  W,D,  s  +  w, 


Wl  =  diag{ej2nfh ,  eilnflh ,  ...  ,  ej2*ft»} 

A/  ~  ’  ^(^2)  ’  •••  ’  T  / (G )] 

W/  =  K(h),  w,(t2),  ...  ,  w;(^)]r 


where  srl  is  N  samples  of  the  received  signal  at  /th  sensor,  s  is  N  samples  of  the  transmitted 
signal,  f  is  the  Doppler  shift  and  Dt  is  the  time  sample  shift  operator  by  n,  =  (rl  / resamples. 
We  can  write  D  -  D"'  where  D  is  an  Ax  A  pennutation  matrix  defined  as  \D\j  =1  if  i  =  j  +  1  , 
[J>W-  !  =  1  and  \D\j  =  0  otherwise. 

According  to  [9]  and  [19],  the  estimated  transmitter’s  position  in  TDOA/FDOA-based  one- 
stage  method  is  found  as  follows.  Let  {pj^be  grid  points  of  possible  emitter  locations.  For 
each  grid  point  form  the  matrix 


DhWhs 

fr\  Ai  ’ 


D  HWHs 


DLHWLHsrL 


5 


where  the  dalay  and  Doppler  operators  correspond  to  the  path  between  the  grid  point  and  the 
respective  sensors  [9],  For  each  grid  point,  then  form  the  LxL  matrix  Qp  =  Vp  Vp  and  find  its 

largest  eigenvalue.  The  location  estimate  is  the  grid  point  that  maximizes  the  largest  eigenvalue, 
that  is 

p  =  arg  max  J/lmax  ( Qp  ) } .  (3) 


Since  all  of  this  process  should  be  done  for  each  one  of  the  grid  points,  it  is  clear  that  a  large 
amount  of  computation  must  be  done  to  find  the  location  and  all  of  it  is  done  at  the  common 
sensor. 


Approved  for  Public  Release;  Distribution  Unlimited. 
5 


It  is  interesting  to  mention  that  in  TDOA/FDOA  based  one-stage  method,  the  (ij) th  element 
of  the  matrix  Q  is  the  value  of  CAF  between  the  signals  received  by  sensors  i  and  sensor  j  [9], 
[19]  and  that  is  why  in  [19],  this  matrix  is  called  the  Cross  Ambiguity  Matrix  (CAM): 


[Qlj  =  W 


y\ = s, 


iHWiD,DjHWjHs. 


=  \CAFy\ 


(T,a)  ■ 


(4) 


where  r  and  co  are  the  corresponding  TDOA  and  FDOA  between  sensors  i  and  j  and  emitter 
located  at  a  specific  emitter  position.  Note  that  the  diagonal  elements  in  CAM  are  the  Auto 
Ambiguity  Function  of  the  received  signals  at  TDOA  =  0  and  FDOA  =  0  which  is  equal  to  the 
energy  of  the  received  signals. 

In  the  following,  we  develop  some  methods  to  increase  the  flexibility  and  feasibility  and 
improve  the  performance  of  one-stage  geolocation  method. 


4,1.2  Approximated  DPD 

As  mentioned  above,  in  one-stage  geolocation  method,  we  need  all  the  received  signals 
together  to  start  the  location  estimation  process.  To  achieve  that,  all  sensors  need  to  transmit 
their  received  signals  to  a  common  point.  The  common  site  performs  a  huge  mathematical 
computation  on  the  raw  received  signals  to  fonn  a  series  of  matrices  used  in  location  estimation 
[7], [8],  [9].  In  this  section  we  develop  a  method  to  distribute  and  reduce  the  amount  computation 
based  on  the  eigenvalue  approximation. 

Definition  (Gershgorin’s  disc)[29]:  Assume  that  A  is  an  n  x  n  complex  valued  matrix 
with  entries  atj  and  Pt  =  y  is  the  summation  of  the  absolute  values  of  all  non-diagonal 

elements  of  the  /'th  row.  Then,  the  set  Dj  =  j  zeC:  \z  -  au\<  P\  is  called  the  /'th  Gershgorin’s  disc 

of  A ,  where  C  is  the  set  of  complex  numbers.  This  disc  contains  the  interior  and  boundary 
points  of  a  circle  with  radius  of  Pt  and  centered  at  au  in  complex  plane. 

Theorem  1  (Gershgorin’s  Theorem)  [29]:  Every  eigenvalue  of  matrix  A  =  \^aij .J  e  C'mn 

lies  within  at  least  one  of  the  Gershgorin’s  discs.  In  other  words,  every  eigenvalue  A  of  matrix  A 
satisfies: 


V^,3  i  \^-ah\<P, 

*?=Xkl- 


(5) 


Theorem  2  [29]:  Assume  that  A  is  an  n  x  n  complex  valued  matrix  with  entries  atJ , 

R .  =  V  n  is  the  summation  of  the  absolute  values  of  all  elements  in  the  /th  row  and  T.  =  V  \ar\ 

j  i 

is  the  summation  of  the  absolute  values  of  all  elements  in  the  /h  column.  Let  R  =  max  Ri  and 

i 

T  =  max  Tj .  Then,  the  absolute  value  of  each  eigenvalue  A  of  matrix  A  satisfies: 


V/l,  |T|  <  min(i?,r) .  (6) 
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As  mentioned  above,  in  TDOA/FDOA  based  one-stage  method,  the  (7,/)th  element  of  the 
CAM  or  Q  in  (3)  is  the  value  of  Cross  Ambiguity  Function  between  the  signals  received  by 
sensors  i  and  sensor  j  ([9], [19])  and  we  name  it  by  CAFj.  The  emitter  location  is  estimated  by 

computing  the  maximum  eigenvalues  of  the  CAM  (or  Q)  in  each  grid  point.  Since  the  CAM  is 
Hennitian  and  positive  definite,  the  eigenvalues  of  CAM  are  real  and  positive.  Moreover,  since 
CAM  is  a  Hennitian  matrix,  we  have,  R  =  T  in  Theorem2  and  consequently,  min(/C  T)  =  R . 
Thus,  for  the  CAM,  the  inequality  in  (6)  can  be  replaced  by: 

\/A,  A  <  max( CAFf) 

i 

CAF,  =  Y\CAFU\  (7) 

j 

4ax  =  ma  x{CAFt) 

l 

where  Amax  is  the  upper  bound  on  eigenvalues  of  the  CAM. 

Suppose  that  we  have  L  receiving  sensors  and  each  one  of  them  broadcasts  its  received  signal 
to  all  other  sensors  in  the  sensor  network.  Then,  each  sensor  i  is  able  to  compute  all  CA  A.  for 

L 

j=l...  L  and  consequently,  it  is  able  to  compute  CAFj  =  .  Now,  if  we  approximate  the 

j= i 

largest  eigenvalue  of  CAM  by  the  upper  bound  on  the  eigenvalues  ( Amax  «  Amax ),  then  the 
location  estimation  will  be  detennined  by  the  point  having  the  largest  Amax . 

Here  is  the  scenario: 

1-  Each  sensor  broadcasts  its  received  signal  in  the  sensor  network. 

2-  Each  sensor  i  computes  CAFF  s  in  TDOA/FDOA  plane  and  then  maps  them  from 

TDOA-FDOA  plane  to  X-Y  (emitter  position)  plane.  The  mapping  will  be  done  very 
easily  knowing  the  position  and  velocity  of  the  sensors  and  also  the  grid  point 
position. 

L 

3-  Each  sensor  i  computes  CAFi  =  ^  |C4A.  |  by  adding  up  the  C/1  A  ’s  and  then  finds  the 

7=1 

peak  of  CAFj  (named  CAF  peak  )  and  its  location  (xj  peak ,  yt  peak )  and  then  transfers  the 

three  numbers  xi  peak ,  v,  k  and  CAFt  peak  to  a  common  sensor  (or  to  all  other  sensors 

since  there  are  just  three  numbers  and  there  is  no  communication  load  to  transfer 
them).  Note  that  this  step  is  motivated  by  Gershgorin’s  Theorem. 

4-  According  to  (3)  and  (7),  the  emitter  location  estimated  is  taken  as  the  (xi  peak ,  yt  peak ) 

corresponding  to  the  largest  CAFi  peak  over  all 

Note  that  in  the  original  DPD  method,  we  need  to  re-compute  and  form  the  matrix  CAM  (or 
Q)  for  each  grid  point  and  find  the  largest  eigenvalue  of  that  matrix  each  time,  which  leads  to  a 
huge  amount  of  computation  especially  when  the  number  of  receiving  sensors  gets  larger. 
Moreover,  all  of  these  computations  would  be  done  at  one  single  point.  But,  the  method  outlined 
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above  does  not  need  to  fonn  the  matrix  CAM  (or  Q)  at  all  nor  does  it  need  to  do  computationally 
expensive  computations  of  the  largest  eigenvalue  each  time.  Thus,  in  the  new  method,  not  only 
has  the  costly  eigenvalue  computation  been  removed,  but  also  the  process  is  distributed  among 
all  receiving  sensors.  When  implementing  DPD  in  a  scenario  where  each  sensor  has  limited 
computational  abilities  it  is  desirable  to  minimize  the  amount  of  computation  done  by  each 
sensor  rather  than  minimize  the  total  computational  complexity.  To  compare  the  computational 
load  suppose  that  there  are  L  sensors  trying  to  estimate  the  emitter  location  in  an  NxN  grid 
plane.  In  the  original  DPD  method,  the  common  sensor  needs  to  compute  the  CAM  for  each  grid 
point.  Since  CAM  is  a  Hennitian  L  xL  matrix  formed  by  CAFs,  the  common  sensor  just  needs  to 

find  all  the  entries  on  and  above  the  main  diagonal.  This  is  equivalent  to  computing  y  (L  - 1) 

CAFs  (as  non-diagonal  elements)  and  L  signal  energies  (as  diagonal  elements).  Moreover,  the 
common  sensor  needs  to  calculate  the  largest  eigenvalue  of  the  matrix  for  each  grid  point  ( N 2 
times).  On  the  other  hand,  in  the  suggested  method,  each  sensor  just  needs  to  find  (L  - 1)  CAFs 
and  one  signal  energy.  In  addition,  they  don’t  need  to  fonn  the  matrix  CAM  and  find  its 

eigenvalues.  Thus,  rather  than  having  one  sensor  compute  —  (Z-l)  CAFs  as  in  the  original 

DPD,  in  the  method  propose  here  each  sensor  computes  only  (L  -  1)  CAFs;  furthermore,  each 
sensor  performs  a  simple  Gershgorin  estimation  rather  than  a  complex  eigenvalue  computation. 

It  is  worth  saying  that  in  the  proposed  method,  if  we  lose  anyone  of  the  sensors  or  even  if  we 
lose  a  couple  of  them,  it  may  reduce  the  accuracy  of  estimation  because  of  missing  some  data 
but,  the  rest  of  the  receivers  can  continue  the  estimation  process  with  no  interruption. 

In  the  proposed  method,  if  we  ignore  the  maximum  operator  term  in  (7)  and  just  take 

Amax  =  CAFt  for  only  one  arbitrary  sensor  i  ,  then  the  results  will  be  equivalent  to  a  method 
named  CAF-MAP  in  [6]  which  has  less  quality  compared  to  DPD  and  Approximated  DPD.  In 
[19],  we  can  also  see  another  approach  named  as  pair-wise  maximum  CAF  detector  that  is  based 
on  comparing  the  value  of  max  max \CAF:j  with  a  threshold  yCAF  for  only  one  arbitrary  i  as 

reference  sensor  (  t,oj  are  TDOA  and  FDOA).  The  results  in  [19]  showed  that  this  method  also 
has  much  lower  quality  in  detection  compared  to  the  GLRT  detector  based  on  largest  eigenvalue; 
no  results  were  provided  in  [19]  on  the  location  accuracy  of  the  pair-wise  maximum  CAF 
method. 

The  simulation  results  for  many  different  cases  show  that  the  eigenvalue  upper  bound  is  very 
close  to  the  true  largest  eigenvalue.  However,  this  approximation  lowers  the  quality  of  the 
estimation  slightly.  We  examined  the  effect  of  the  proposed  approximation  on  the  estimation 
accuracy  using  Monte-Carlo  computer  simulations  (with  500  runs  each  time).  In  this  simulation, 
a  set  of  8  moving  sensors  and  one  stationary  emitter  are  placed  in  a  configuration  as  shown  in 
Figure  1.  There  exists  a  cross  ambiguity  function  for  each  two  of  the  sensors.  The  sampling 
frequency  is  80  kHz  and  the  number  of  samples  is  equal  to  4096.  Figure  2  shows  the  effect  of 
eigenvalue  approximation  on  RMS  error  of  emitter  location  estimation  for  X  and  Y  dimensions. 
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Figure  1 :  Placement  of  the  sensors  and  the  emitter  position  used  for  simulation. 


(a)  (b) 


Figure  2:  RMS  errors  for  X  and  Y  versus  SNR. 


In  the  suggested  method,  if  we  ignore  the  maximum  operator  term  in  (7)  and  just  take 
Amax  =  CAF.  for  only  one  arbitrary  sensor  i  ,  then  the  results  will  be  equivalent  to  a  method 

named  CAF-MAP  in  [6]  which  has  less  quality  compared  to  DPD  and  Approximated  DPD.  In 
[19],  we  can  also  see  another  approach  named  as  pair-wise  maximum  CAF  detector  that  is  based 
on  comparing  the  value  of  max  max  CAR.  with  a  threshold  yCAF  for  only  one  arbitrary  i  as 

j  T,0)  I  J  I 

reference  sensor  (  t,co  are  TDOA  and  FDOA).  The  results  in  [19]  showed  that  this  method  also 
has  much  less  quality  in  detection  compared  to  the  GLRT  detector  based  on  largest  eigenvalue. 
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4.1.3  Decentralized  DPD 

In  this  section,  we  develop  another  method  to  implement  DPD  with  the  goal  of  reducing  the 
computation  load  on  the  center  point  by  distributing  the  mathematical  computations  among  all 
receivers  and  using  data  compression  to  reduce  the  amount  of  data  transmission.  In  this  method, 
we  don’t  apply  any  approximation  more  than  data  compression.  Thus,  the  accuracy  of  this 
method  is  as  same  as  original  DPD  for  high  enough  bit  rates. 

4.1.3. 1  Distributed  Computation 

As  mentioned  above,  in  DPD  method  a  common  site  perfonns  a  huge  mathematical 
computation  on  the  raw  signals  received  from  all  other  sensors  to  form  a  series  of  matrices  used 
in  location  estimation.  If  each  sensor  or  pair  of  sensors  can  do  some  pre-processing  on  its  own 
data  before  transmitting  to  the  common  site,  it  can  help  to  do  the  estimation  in  shorter  time  and 
with  less  processing  load  on  the  common  site.  It  also  helps  to  reduce  the  sensitivity  of  the 
common  site’s  role  and  the  dependence  of  the  whole  process  to  one  single  point  in  the  sensor 
network. 

As  mentioned  above,  in  TDOA/FDOA  based  one-stage  method,  the  (i,j) th  element  of  the 
cross  ambiguity  matrix  (CAM  in  [20]  or  Q  in  [9]  and  (3))  at  grid  point  p  is  the  value  of  Cross 
Ambiguity  Function  (CAF)  between  the  signals  received  by  sensors  i  and  sensor  j  at 
TDOA/FDOA  point  corresponding  to  the  grid  point  p.  Thus,  sensor  i  and  sensor  j  can  contribute 
to  compute  the  (ij)th  element  of  CAM  (and  (j,i) th  element  as  well  since  CAM  is  a  Hennitian 
matrix).  In  the  other  word,  each  pair  of  sensors  i  and  j  can  share  their  received  signals  together  to 
compute  the  CA  Fy  in  TDOA/FDOA  domain.  Then  the  computed  CAFs  can  be  transmitted  to  the 
common  site.  The  common  site  only  maps  the  received  CAFs  to  X-Y  plane  and  uses  them  to 
form  the  CAM  matrix. 

Figure  3  shows  a  simple  case  with  three  receiving  sensors.  In  Figure  3(a),  we  can  see  the 
sensors  sharing  their  received  signals.  Sensor  1  sends  its  received  signal  to  sensor  2,  sensor  2 
sends  its  received  signal  to  sensor  3  and  sensor  3  sends  its  received  signal  to  sensor  1.  Then, 
sensor  1,  sensor  2  and  sensor  3  are  able  to  compute  CAFn,  CAFn,  CAFn  respectively.  The 
tenns  Au,  A 22,  A 23  are  the  Auto  Ambiguity  Function  of  the  received  signals  at  the  origin  of 
TDOA/FDOA  plane  that  are  equal  to  the  energy  of  the  received  signals  at  each  receiver.  In 
Figure  3(b),  we  can  see  the  receivers  transmitting  the  computed  CAFs  and  Energies  to  a  common 
site  (that  obviously  can  be  any  one  of  the  three  sensors)  to  fonn  the  matrix  CAM.  Note  that  each 
one  of  the  CAF13,  CAFn  and  CAF 23  is  a  complex-valued  matrix  used  to  fill  out  the  off-diagonal 
elements  of  the  CAM  and  An,  A 22,  A 33  are  just  three  real  numbers  sitting  on  the  diagonal 
elements  of  the  CAM. 
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Figure  3:  Decentralized  DPD  for  three  receivers. 


In  this  method,  all  sensors  contribute  in  the  process  of  estimation  and  the  computation  load  is 
divided  among  all  members  of  the  sensor  network.  However,  this  method  does  not  help  to 
eliminate  the  bottle-neck  we  have  in  data  transmission  since  we  need  to  transmit  all  CAFs  to  the 
common  point.  In  next  sub-sections,  we  develop  some  data  compression  ideas  to  compress  the 
CAF  before  transmission  to  the  common  site  to  reduce  the  amount  of  data  transmission.  We 
exploit  some  of  the  CAF  properties  to  eliminate  the  redundancy  and  obtain  a  larger  compression 
rate. 


4.1.3.2  Exploiting  CAF  Properties  for  Data  Compression 

To  reduce  the  amount  of  data  transmission  and  channel  bandwidth  in  our  sensor  network  and 
obtain  a  better  perfonnance,  we  apply  some  data  compression  methods  in  different  levels  of  data 
transmission.  In  [4]  and  [5],  Fowler  and  Chen  have  derived  some  fisher  information  based 
methods  to  compress  the  received  signal  in  each  pair  of  sensors.  Here,  we  develop  and  apply 
some  proper  compression  methods  on  computed  CAFs  before  transferring  it  to  the  center  site. 

CAF  is  a  two-dimensional  complex  valued  function.  Thus,  we  can  consider  the  CAF  to  be  an 
image  (albeit  a  complex  valued  image)  and  apply  image  compression  methods  to  it.  For  example, 
Embedded  Zerotree  Wavelet  (EZW)  [10]  can  be  used  to  compress  the  CAF.  One  of  the  most 
important  reasons  that  encourage  us  to  use  EZW  method  is  that  EZW  is  an  embedded  algorithm. 
It  attempts  to  provide  a  sequence  of  bits  that  if  truncated  anywhere  gives  the  best  distortion  for 
that  rate.  The  EZW  compression  causes  two  kinds  of  distortion  on  the  data.  The  first  distortion 
refers  to  the  effect  of  data  quantization  which  can  be  modeled  by  an  additive  noise.  The  second 
distortion  is  the  result  of  throwing  away  the  insignificant  coefficients  of  the  wavelet  transform 
which  can  be  roughly  modeled  as  passing  the  data  through  a  lowpass  filter.  Note  that  a  typical 
CAF  contains  a  large  main  lobe  and  (usually)  various  small  side  lobes.  An  important  aspect  for 
compression  of  the  CAF  is  that  it  is  a  relatively  slowly  changing  function.  The  fast  changing 
parts  (which  are  equivalent  to  very  high  frequency  points)  come  from  the  effect  on  the  CAF  of 
the  additive  noise  of  received  signals.  Thus,  viewed  as  an  image,  it  seems  that  the  important  part 
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to  be  retained  is  a  spatially  low-pass  type  signal  that  should  show  up  in  the  medium  and  low 
frequency  parts  of  the  wavelet  transform  used  in  EZW.  Because  most  of  the  data  concentrated  in 
medium  and  low  frequency  parts,  it  helps  to  obtain  lots  of  zero  tree  roots  and  that  is  another 
encouraging  reason  to  use  EZW  method.  Because  the  high  frequency  parts  of  the  wavelet 
transform  will  contain  mostly  noise  and  will  be  discarded  (except  at  the  very  highest  data  rate, 
when  little  or  no  truncation  is  done),  the  EZW  algorithm  will  also  perform  a  denoising  operation. 

There  are  some  other  special  properties  with  CAF  that  can  be  exploited  to  get  better 
compression  performance.  One  of  the  most  important  properties  of  CAF  is  the  symmetry.  In  the 
Appendix  we  have  proved  that  CAF  is  symmetric  around  its  magnitude  peak: 

\CAFy(-T  +  Tp,-(o+o)p)\  =  {CAF^t +  tp,co  +  cop)\  (8) 

where  CAFtJ  is  the  cross  ambiguity  function  between  the  signals  received  by  sensor  i  and  sensor  j 
and  ( jp ,  o>p )  is  the  peak  of  CAF  magnitude.  This  result  provides  a  kind  of  symmetry  of  the  CAF 
around  the  point  (rp,  cop)  or  the  peak  of  CAF  magnitude  which  can  be  exploited  for  data 
compression.  In  practice,  the  received  signals  received  signals  are  the  delayed  and  Doppler- 
shifted  version  of  transmitted  signal  plus  noise.  This  noise  perturbs  the  CAF  a  little  bit  from  the 
perfect  symmetry.  Thus,  we  rewrite  (8)  as, 

\CAF.j(-t  +  tp,-g)  +  (qp)\  =  \  CAF.j(T  +  Tp,eo  +  a>p)\  + E  (9) 

where  E  can  be  the  error  from  perfect  symmetry  which  is  a  negligible  value.  Thus,  using  the 
symmetry  property,  it  is  possible  to  extract  the  entire  CAF  magnitude  by  transmission  of  only 
half  of  the  CAF  magnitude  plus  the  small  residual  amount  of  E. 

We  examined  the  performance  of  the  proposed  method  using  Monte  Carlo  computer 
simulations.  In  this  simulation,  a  set  of  4  moving  sensors  and  one  stationary  emitter  are  placed  in 
a  configuration  as  shown  in  Figure  4.  There  exists  a  cross  ambiguity  function  for  each  two  of  the 
sensors  that  should  be  compressed  and  transmitted  to  a  common  site  to  do  the  location 
estimation.  The  signals  are  Binary  Phase  Shift  Keying  (BPSK)  signals,  the  sampling  frequency  = 
400  kHz,  SNR  =  -10  dB  and  the  number  of  samples  is  equal  to  65536. Two  different  compression 
methods  have  been  examined  in  this  simulation.  In  the  first  method,  we  just  applied  the  EZW 
algorithm  to  compress  the  CAF  (labeled  “simple  compression”  in  the  figures).  In  the  second 
method,  we  applied  the  EZW  algorithm  to  compress  the  CAF  and  we  used  the  symmetry 
property  to  reduce  the  amount  of  transmitted  data  (labeled  “symmetric  compression”  in  the 
figures).  The  effect  of  data  compression  on  RMS  error  of  emitter  location  estimation  for  X  and  Y 
dimensions  is  illustrated  in  Figure  5  (a)  and  (b),  respectively.  Obviously,  the  RMS  error  will 
decrease  by  increasing  the  bit  rate.  Comparing  the  two  curves  in  each  plot  shows  that  the 
symmetric  compression  method  gives  us  more  accurate  results  for  the  same  bit  rates. 
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Figure  4:  Placement  of  the  sensors  and  the  emitter  position  used  for  simulation. 


Figure  5:  RMS  errors  for  X  and  Y  versus  bits/element. 


4.1.3.3  SVD-based  Data  Compression  for  CAF 

The  singular  value  decomposition  (SVD)  is  an  important  tool  with  many  useful  signal 
processing  applications.  For  a  complex  valued  M  xN  matrix  X,  the  SVD  representation  will  be 

X  =  UEV^  =  Vcruvf 

l  l  l 

1=1 
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where  U  is  an  M  x  M  unitary  matrix  consisting  of  M  left  singular  vectors  as  its  columns,  V  is 
an  N  xN  unitary  matrix  consisting  of  N right  singular  vectors  as  its  columns  and  £  is  a  pseudo¬ 
diagonal  M  xN  matrix  with  nonnegative  real  singular  values  (07)  on  the  main  diagonal  ordered 
such  that  cr  >  cr.+l ,  r  is  the  number  of  non-zero  singular  values,  u,.  is  the  7th  left  singular  vector 

and  vf  is  the  Hennitian  transpose  of  the  7th  right  singular  vector.  By  truncating  the  above 
summation  to  k  <  r  terms,  we  get  a  rank-A:  matrix 


x,=UAv"  =  Z<wf 


that  approximates  X  better  than  any  other  rank-A:  matrix  in  the  least  square  error  sense  [23],  [24]. 
This  is  the  main  idea  of  SVD  data  compression. 

The  complex-valued  M  xN  matrix  X  contains MN complex  values  or  equivalently  2MN real 
values.  However,  the  truncated-SVD  version  Xk  uses  only  kM  complex  values  to  represent 
matrix  Ut ,  kN  complex  values  for  matrix  Yk  ,  and  k  real  values  to  represent  the  singular  values. 
Thus,  in  approximation  X  by  X*,  the  compression  ratio  is: 


CR 


2MN 

2  kM  +  2  kN  +  k  ' 


For  example,  for  a  128  x32  matrix  truncated  for  k=  1  the  compression  ratio  is  25 : 1 . 

As  mentioned  in  [25],  the  singular  value  decomposition  of  an  image  is  conceptually  similar 
to  its  Karhunen-Loeve  decomposition  but  in  a  different  manner.  The  first  difference  is  that 
Karhunen-Loeve  decomposition  basis  are  determined  by  the  covariance  matrix  of  the  random 
process  that  generates  the  image  but,  SVD  is  defined  on  the  image  itself.  The  second  difference 
is  that  if  both  representations  are  truncated  for  the  purpose  of  data  compression,  SVD  is  the  best 
approximation  in  least  square  error  sense,  while  Karhunen-Loeve  is  the  best  approximation  in 
mean  square  error  sense. 

CAF  usually  contains  a  big  main  lobe  and  several  small  side  lobes  that  if  we  slice  each  of 
them  up  at  different  points,  we  will  always  get  a  curve  with  a  similar  shape.  It  has  been  shown 
that  for  a  time-frequency  localization  operator  there  are  several  large  singular  values  at  the 
beginning,  followed  by  a  sharp  plunge  in  the  values,  with  a  final  asymptotic  decay  to  zero  [26]. 
Since  the  cross  Ambiguity  function  is  considered  to  be  a  member  of  Cohen’s  class  of  time- 
frequency  representations  [27],  these  properties  imply  that  CAF  is  very  close  to  a  low  rank 
matrix.  Thus,  most  of  the  data  is  concentrated  in  the  first  few  singular  vectors  and  values. 

In  reality,  the  received  signals  are  noisy.  The  effect  of  the  noise  on  the  singular  values  is 
spread  among  all  the  singular  values  but,  as  mentioned  before,  most  of  the  data  is  concentrated  in 
the  first  few  singular  vectors  and  values.  Thus,  by  SVD  truncation  we  reduce  the  amount  of  noise 
and  equivalently  we  increase  the  SNR  [28].  The  singular  values  of  a  sample  128  x32  CAF  are 
illustrated  in  Figure  6  for  two  cases:  (a)  noiseless  signals  and  (b)  noisy  signals.  As  we  can  see, 
there  are  only  3  to  5  significant  singular  values  in  Figure  6(a)  showing  that  the  CAF  is  very  close 
to  a  low  rank  matrix.  But,  Figure  6(b)  shows  that  in  the  noisy  case  the  number  of  significant 
singular  values  increases  to  12.  Therefore,  it  is  clear  that  the  signal  to  noise  ratio  can  increase  by 
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applying  SVD  data  compression  and  retaining  the  first  few  singular  values  and  discarding  the 
rest. 


(a)  (b) 


Figure  6:  Singular  Values  of  CAF:  (a)  Noiseless  signal,  (b)  Noisy  signal 


We  have  applied  an  SVD-based  data  compression  on  CAF.  The  results  show  that  SVD 
approach  is  a  beneficial  method  for  CAF  data  compression  and  also  it  is  a  strong  tool  for  CAF 
denoising.  Simulation  results  show  that  by  applying  SVD  Data  Compression  it  is  possible  to 
perfonn  accurate  location  estimation  in  spite  of  the  fact  that  we  transmit  fewer  bits.  Also  for 
smaller  compression  ratio  (higher  bit  rates),  we  even  achieve  an  improvement  in  performance  of 
location  estimation  compared  to  the  case  that  we  do  not  compress  the  data  at  all  and  that  is 
because  of  the  de-noising  effect  of  the  SVD. 

We  examined  the  performance  of  the  proposed  method  using  Monte  Carlo  computer 
simulations  (with  500  runs  each  time).  In  this  simulation,  a  set  of  4  moving  sensors  and  one 
stationary  emitter  are  placed  in  a  configuration  as  shown  in  Figure  7.  There  exists  a  cross 
ambiguity  function  for  each  two  of  the  sensors  that  should  be  compressed  and  transmitted  to  a 
common  site  to  do  the  location  estimation.  In  this  simulation,  the  signals  are  BPSK,  the  sampling 
frequency  =  20  kHz  and  the  number  of  samples  is  equal  to  4096.  Figure  8  shows  the  effect  of 
data  compression  on  RMS  error.  The  four  curves  compare  the  cases  (i)  without  compression,  (ii) 
SVD-based  compression  with  compression  ratio  of  25:1,  (iii)  SVD-based  compression  with 
compression  ratio  of  8:1,  and  (iv)  SVD-based  compression  with  compression  ratio  of  5:1.  As  we 
can  see,  even  for  high  compression  ratio  of  25:1,  the  estimation  accuracy  is  pretty  close  to  the 
case  without  compression.  Surprisingly,  the  case  with  the  compression  ratio  of  5 : 1  (and  even  the 
case  with  the  compression  ratio  of  8:1  in  some  points)  yields  more  accurate  results  than  without 
compression  case.  This  improvement  is  obtained  because  of  the  de-noising  property  of  SVD- 


Approved  for  Public  Release;  Distribution  Unlimited. 
15 


based  data  compression. 
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Figure  7:  Placement  of  the  sensors  and  the  emitter  position  used  for  simulation. 


(a)  (b) 


Figure  8:  Location  error  using  SVD-based  CAF  compression  for  various  compression  ratios. 
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4.1.4  Semi-Decentralized  DPD 

Comparing  the  classic  method  to  DPD  method,  we  can  see  that  in  classic  method  we  only 
transfer  the  TDOA/FDOA  values  to  the  center  site  that  leads  to  a  very  low  data  transmission 
load;  but  in  DPD  we  transfer  the  whole  received  signals  (or  the  whole  CAFs  as  suggested  in 
previous  section)  to  the  center  site  that  ends  up  with  more  accurate  results  in  localization.  In  this 
section,  we  develop  a  method  taking  the  advantages  of  both  classic  and  DPD  scenarios.  In  fact, 
we  use  the  same  idea  of  Decentralized  DPD  exploiting  the  relationship  between  different  CAFs 
to  reduce  the  amount  of  data  transmission  in  the  sensor  network 

Suppose  that  two  sensors  R1  and  R2  receive  the  LPE  signals  u(t)  and  v(t),  respectively. 
Under  the  so-called  narrowband  approximation  and  assuming  that  we  can  estimate  the  antenna 
and  electronic  devices  attenuations,  the  lowpass  equivalent  (LPE)  model  of  the  received  signal 
for  noise-free  case  will  be  [11]: 


u(t)  =  ae  M's(t  -  r, ) 
v(t)  =  ae~j0}2‘s(t  ~t2) 

where  r,  and  z2  are  the  time  delays  and  o\  and  ox  are  the  Doppler  shifts  for  the  first  and 
second  received  signals,  a  =  e~JCOcT'  ,  (3  =  A7*cT2  and  coc  is  the  carrier  frequency  in  rad/sec  [11]. 
Now,  we  can  write  one  of  the  received  signals  in  terms  of  the  other  one: 

v(0  =  u(t  +  rn)eirajneM2!eMT'2  (10) 

where  r12  =  -  r2)  is  the  TDOA  and  con  =  (cox  -  co2 )  is  the  FDOA.  Then,  the  CAF  is: 

oo 

CAF12(t,co)  =  |  u{t)v*{t-  r)eJwtdt 

-co 

oo 

=  eM2r-MTl2~MTl2  J  u(t)u(t  -  (-r12  +  rjy^-^dt  (11) 

-00 

=  eM2T-Ml'2~Mr'2  An(r-Tl2,(Q-con) 


where  the  second  line  follows  from  (10)  and  the  third  line  follows  after  defining 


An{r,(o)  =  ^u{t)u  {t-r)ejMdt  , 

-oo 

which  is  the  Auto  Ambiguity  function  (AAF).  With  a  similar  approach  for  three  sensors  and 
applying  (10)  for  the  signals  received  by  first  and  second  receivers,  we  will  have: 

CAF  12  ( r ,  (6)  =  g^(-®-®i2+®c+®i  )ra  CAFn(z  +  r12 ,  a  +  con )  (12) 

where  CAFmn  is  the  cross  ambiguity  function  between  signals  received  by  sensor  m  and  sensor  n. 
In  equation  (12),  having  (r12,o12)  and  one  sample  point  of  CAF2i  and  CAF i3  it  is  possible  to 
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compute  cox .  In  the  noisy  case,  the  point  (rl2,ry|2)  can  be  estimated  as  the  position  of  peak  of  the 

CAF  12  magnitude  in  the  (TDOA-FDOA)  plane  and  that  is  why  we  named  this  method  as  Semi 
Optimal  Decentralized  DPD.  Also,  note  that  this  estimation  is  applicable  as  long  as  we  are  able 
to  find  the  peak  of  the  CAF  12  magnitude.  Thus,  having  (rn,a>xl),  CAF  13  ,  it  is  possible  to 

reconstruct  the  whole  CAF23  and  we  don’t  need  to  transmit  the  CAF23  to  the  center  site  anymore. 
Similarly,  in  a  sensor  network  when  L  sensors  receive  a  signal  from  a  single  emitter,  we  can 
show  that: 


CAF„m{T,(0)  =  -’^’A’-CAFJt  +  ffl,)  (13) 

where  t  and  a  are  the  TDOA  and  FDOA  between  signals  received  by  sensor  p  and  sensor 
m.  As  a  special  case  when  p  =  1,  we  have: 

=  +  (14) 

Exploiting  (14),  it  is  possible  to  construct  the  Cross  Ambiguity  Matrix  (CAM  in  [19]  or  Q  in 
[9])  by  transferring  only  the  first  element  in  each  column  of  the  matrix  (in  other  word,  by 
transferring  the  first  row  of  the  matrix)  plus  the  corresponding  TDOA/FDOAs  and  it  leads  to  a 
large  reduction  in  computation  load  and  amount  of  data  transmission. 

Figure  9  illustrates  a  simple  case  of  4  receiving  sensors.  As  we  see  in  Figure  9(a),  sensorl  (it 
can  be  obviously  each  of  the  sensors)  propagates  its  received  signal  into  the  network.  Each  one 
of  the  other  sensors  can  compute  the  CAF  between  sensorl  and  itself.  Then,  each  sensor  can 
transmit  the  computed  CAF  to  a  common  site  (which  can  be  one  of  the  sensors).  In  Figure  9(b), 
we  assumed  sensor2  as  the  common  point.  After  receiving  CAF  13  and  CAF  14  ,  sensor2  is  able  to 
compute  CAF23,  having  CAFn  (computed  by  itself)  and  CAF  13  (received  from  sensor3)  applying 
(14).  It  is  also  able  to  compute  CAF2 4  and  CAF34  having  CAFn,  CAF  13  and  CAF14  easily.  Now, 
sensor2  has  everything  to  form  the  matrix  CAM  and  do  the  location  estimation.  Note  that 
according  to  (14),  all  of  the  CAF  calculations  are  nothing  more  than  delay,  Doppler  and  phase 
shifts.  Thus,  the  sensor2  does  not  suffer  from  extra  CAF  calculations. 
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We  examined  the  proposed  method  using  Monte-Carlo  computer  simulations  (with  500  runs 
each  time).  In  this  simulation,  a  set  of  3  moving  sensors  and  one  stationary  emitter  are  placed  in 
a  configuration  as  shown  in  Figure  10.  In  this  simulation,  we  compute  CAFn  and  CA F\ 3  directly 
using  stein  method  [2]  and  then  compute  the  CAF23  using  (12).  Then  having  all  of  the  elements, 
we  formed  matrix  CAM  and  apply  the  DPD  method.  The  sampling  frequency  is  20  kHz  and  the 
number  of  samples  is  equal  to  4096.  Figure  1 1  shows  the  RMS  error  of  emitter  location 
estimation  for  X  and  Y  dimensions  compared  to  original  DPD  and  Classic  methods. 
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Figure  10:  Placement  of  the  sensors  and  the  emitter  position  used  for  simulation. 


Semi  Decentralized  DPD  leads  to  a  large  reduction  in  computation  load  compared  to  DPD 
and  large  reduction  in  the  amount  of  data  transmission  compared  to  Decentralized  DPD  with  the 
cost  of  lower  accuracy,  but  as  we  can  see  in  Figure  1 1,  it  is  still  much  more  accurate  compared  to 
classic  TDOA/FDOA  method. 
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Figure  1 1 :  RMS  errors  for  X  and  Y  versus  SNR. 


4.1.5  Summary 

The  one-stage  localization  method  is  a  major  new  development  in  TDOA/FDOA-based 
emitter  location  that  provides  significant  improvement  in  performance  at  low  SNR  levels. 
However,  that  improvement  comes  at  a  cost  of  significantly  more  computational  complexity. 
Worse,  as  DPD  was  proposed,  that  complexity  is  all  concentrated  at  one  computing  node,  which 
is  different  from  the  classical  method  where  the  computations  are  distributed  evenly  among  the 
sensors.  Furthermore,  unlike  the  classical  method,  the  location  processing  is  highly  complex 
(requiring  the  computation  of  eigenvalues  for  each  grid  point).  The  large  amount  of 
mathematical  computations  that  should  be  done  by  a  single  center  point  may  make  some 
restrictions  in  processing.  Moreover,  a  complete  dependence  of  the  whole  process  (in  both  data 
collection  and  data  processing)  on  only  one  single  point  is  also  a  problem  that  may  reduce  the 
reliability  of  the  system.  In  this  paper,  we  developed  several  methods  to  reduce  the  load  of  data 
computation  and  data  transmission  using  distributed  data  computation  and  processing,  applying 
data  compression  methods,  exploiting  the  CAF  properties,  taking  advantage  of  CAF 
relationships  in  the  sensor  network  and  exploiting  some  kind  of  beneficial  approximations. 

The  Approximated  DPD  method  proposed  here  exploits  the  simplicity  of  Gershgorin’s 
theorem  to  approximately  compute  the  largest  eigenvalue  without  the  high  cost  of  exactly 
computing  it.  This  enables  each  sensor  to  locally  make  its  best  estimate  of  the  location  based  on 
that  data  it  has.  These  locally-generated  estimates  are  then  transmitted  to  a  central  location 
where  a  final  decision  is  made. 

In  Decentralized  DPD  method,  we  applied  the  distributed  computation  idea  to  divide  the 
mathematical  calculation  load  among  all  sensors  of  the  network.  In  this  method,  we  don’t  use 
any  other  approximation  more  than  data  compression  and  the  quality  of  this  method  is  the  same 
as  the  Original  DPD  for  proper  bit  rates. 
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Finally,  in  Semi-Optimal  Decentralized  DPP,  we  used  the  same  idea  of  Decentralized  DPD 
in  addition  to  exploiting  the  relationship  between  different  CAFs  to  eliminate  the  data 
redundancy  and  this  idea  leads  to  a  large  data  transmission  reduction. 

All  of  the  proposed  methods  allows  DPD  to  be  implemented  in  a  decentralized  manner  where 
no  single  sensor  is  required  to  do  an  unfair  share  of  the  computations,  yet  the  perfonnance 
improvement  of  DPD  is  not  sacrificed. 


4.2  Spatial  Sparsity-Based  Approach  to  Emitter  Location 

Passive  emitter  localization  is  a  challenging  issue  in  statistical  signal  processing.  The 
position  can  be  estimated  by  measuring  one  or  more  location-dependent  signal  parameters.  One 
of  the  most  popular  and  common  emitter  location  methods  is  based  on  TDOA  and  FDOA 
estimations.  In  the  classical  approach  to  this  method,  FDOA  and  TDOA  are  estimated  from  the 
cross-correlation  of  the  signals  received  by  several  pairs  of  sensors  [1];  this  is  done  by 
computing  the  cross  ambiguity  function  (CAF)  [2]  and  finding  the  peak  of  its  magnitude  surface. 
Then  these  TDOA/FDOA  estimates  are  used  in  statistical  processing  to  locate  the  emitter  [3]. 

However  the  classic  two-stage  method  is  not  necessarily  optimal  because  in  the  first  stage  of 
these  methods,  the  TDOA  and  FDOA  estimates  are  obtained  by  ignoring  the  fact  that  all 
measurements  should  be  consistent  with  a  single  emitter  location  [9].  In  other  words,  each  stage 
is  itself  optimal  but  the  cascade  of  the  two  stages  is  not  necessarily  optimal. 

In  this  section,  we  exploit  spatial  sparsity  of  the  emitter  on  the  x-y  plane  and  use  convex 
optimization  theory  to  estimate  the  location  of  the  emitter  directly  without  going  through  the 
intennediate  stage  of  TDOA/FDOA  estimation.  It  is  obvious  that  in  emitter  location  problems, 
the  number  of  emitters  is  much  smaller  than  the  number  of  all  grid  points  in  a  fine  grid  on  the  x-y 
plane.  Thus,  by  assigning  a  positive  number  to  each  one  of  the  grid  points  containing  an  emitter 
and  assigning  zero  to  the  rest  of  points,  we  will  have  a  very  sparse  grid  plane  matrix  that  can  be 
refonned  as  a  sparse  vector.  Since  each  element  of  this  vector  corresponds  to  one  grid  point  in 
the  x-y  plane,  we  can  estimate  the  location  of  emitters  by  extracting  the  position  of  non-zero 
elements  of  the  sparsest  vector  that  satisfies  the  TDOA/FDOA  relationship  between  transmitted 
signals  and  received  signals.  In  principle,  sparsity  of  the  grid  vector  can  be  enforced  by 
minimizing  its  /(l-norm  (i.e.,  the  number  of  non-zero  elements  in  the  grid  vector).  However, 
since  the  /0-norm  minimization  is  an  NP-hard  non-convex  optimization  problem,  it  is  very 
common  (e.g  in  compressive  sensing  problems)  to  approximate  it  with  /:  ,  -norm  minimization, 
which  is  a  convex  optimization  problem  and  also  achieves  the  sparse  solution  very  well  [30]. 
Thus,  after  formulating  the  problem  in  tenns  of  the  sparse  grid  vector,  we  can  estimate  this 
vector  by  pushing  sparsity  using  i ,  -norm  minimization  on  the  grid  vector,  subject  to  the 
TDOA/FDOA  relationship  between  the  signals  transmitted  from  the  grid  point  and  the  signals 
received  by  the  sensors. 

In  [31],  the  authors  suggested  a  source  localization  method  based  on  TDOA  in  a  multipath 
channel  exploiting  the  sparsity  of  the  multipath  channel  for  estimation  of  the  line-of-sight 
component.  In  this  method,  the  sensors  don’t  need  to  know  the  information  on  the  specific 
transmitted  symbols  but,  they  require  knowledge  of  the  pulse  shape  of  the  transmitted  signal.  In 
[32],  the  authors  suggested  a  compressive-sensing  based  distributed  target  localization  using 
TDOA.  In  this  method,  each  sensor  approximates  the  transmitted  signal  by  its  own  received 
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signal  mapped  to  each  one  of  the  grid  points.  This  idea  helps  to  reduce  the  amount  of  data 
transmission  in  the  sense  of  distributed  localization  but  it  lowers  the  quality  of  the  estimation 
since  each  sensor  estimates  the  transmitted  signal  just  using  its  own  received  signal.  Also,  each 
sensor  computes  its  own  location  estimation  that  is  not  necessarily  equal  to  other  sensors’ 
estimations.  Weiss  and  Amar  [7],  [8],  [9]  developed  a  single-stage  Least-Squares  method  using 
TDOA  and  FDOA,  named  direct  position  detennination  (DPD).  Vankayalapati  and  Kay  [19]  also 
derived  similar  results  based  on  a  detection  theory  point  of  view;  the  DPD  estimator  was  derived 
as  the  ML  estimator  needed  for  the  generalized  likelihood  ratio  detector.  The  performance  of  the 
DPD  method  is  better  than  the  two-stage  classic  method  (especially  for  low  SNRs).  However,  the 
simulation  results  show  that  DPD  does  not  obtain  accurate  results  in  the  case  of  multipath  or 
multi-emitter  scenarios. 

In  this  paper,  contrary  to  [31]  and  [32],  we  developed  a  method  based  on  both  TDOA  and 
FDOA  to  take  advantage  of  both  delay  and  Doppler  shifts.  Contrary  to  [6],  our  method  does  not 
need  any  knowledge  of  the  transmitted  signal’s  pulse  shape  nor  any  other  a  priori  information. 
Similar  to  [32],  we  exploit  the  grid  point  spatial  sparsity  but,  we  consider  the  transmitted  signal 
as  a  deterministic  unknown  signal  that  will  be  estimated  in  the  sensor  network  using  all  received 
signals.  Similar  to  [19]  and  [9],  we  estimate  the  emitter  location  directly  without  going  through 
the  intennediate  stage  of  TDOA/FDOA  estimation.  However,  the  Monte-Carlo  simulation  results 
show  the  higher  performance  of  the  proposed  method  compared  to  DPD  method  and  classic  two- 
stage  method  especially  in  multipath  scenarios. 

Suppose  that  an  emitter  transmits  a  signal  and  L  sensors  receive  that  signal.  The  complex 
baseband  signal  observed  by  the  /th  sensor  is 

r,(t)  =  als{t-r^)ej2nfl‘  +  w,(t)  (15) 

where  s(t)  is  the  transmitted  signal,  a,  is  the  complex  path  attenuation,  ,  f  is  the  Doppler  shift, 
r!  is  the  signal  delay  and  wl  (t)  is  a  white,  zero  mean,  complex  Gaussian  noise.  Assume  that  each 
sensor  collects  Ns  signal  samples  at  sampling  frequency  Fs  =  1  /  Ts .  Then,  we  have 

rt  =  al  Wl  Dt  s  +  wt  (16) 

with 

s  =  [s(tl),  s(t2),  ...  ,  s(t^)f 
r,=[r,(tx),  r,(t2) ,  ...  ,  r;(^)f 

w)  =  K(h),  w,(t2),  ...  ,  w,(tNs)Y 

W,  =  diag{eilKflh ,  ei2nflh,  ...  ,  eilxf,N°} 

where  r,  is  the  vector  containing  Ns  samples  of  the  received  signal  by  /th  sensor,  s  is  Ns  samples 
of  the  transmitted  signal,  f  is  the  Doppler  shift  and  Z)/  is  the  time  sample  shift  operator  by 
n,  =  (t,  /  7) )  samples.  We  can  write  D,  =  D  where  D  is  an  Ns  x  /Vv  permutation  matrix  defined 
as  [D\j  =1  if  i  =  j  + 1  ,  [Z)]0  v_,  =  1  and  [D]v  =  0  otherwise. 
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Now,  we  assign  a  number  zt  to  each  one  of  the  grid  points  (x,y).  Assume  that  zx  is  one  for 

the  grid  points  containing  an  emitter  and  zero  for  the  rest  of  the  grid  points.  Thus,  the  signal 
vector  received  by  /th  sensor  will  be 


'i=II 


ZX,yai,x,yWl,x,yDl 


,x,ys  +  wn 


(17) 


where  Wlxy  and  Dlx  are  the  Doppler  shift  and  time  sample  shift  operators  w.r.t  sensor  / 

assuming  that  the  emitter  is  located  in  the  grid  point  (x,y)  and  the  summations  are  over  all  grid 
points  in  the  desired  (x,y)  range.  Now,  if  we  reform  all  of  the  grid  points  in  a  column  vector  and 
re-arrange  the  indices,  we  will  have 


ri=TjZnai,nWl,nDl,nS  + 


W, 


n- 1 


(18) 


In  (18),  we  consider  the  transmitted  signal  s  as  a  deterministic  unknown  signal  (a  common 
signal  model  in  localization  problems).  Then,  for  each  grid  point,  we  estimate  the  transmitted 
signal  using  the  Minimum  Variance  Unbiased  estimator  (MVU)  as 


L  1= 1 


(19) 


where  sn  is  the  MVU  estimate  for  the  transmitted  signal  from  grid  point  n. 

We  define  the  matrix  r n  as  the  Doppler  and  delay  operator  w.r.t  all  L  sensors,  assuming  that 
the  received  signal  comes  from  the  grid  point  n  (there  is  an  emitter  at  grid  point  n): 


r  = 


a-,  W,  D1 

2,n  2,n  2  ,n 


a,  W,  D, 

L,n  L,n  L,n  JLN  xn 


Then,  we  can  define  0n,n  e{l,  2,  ...,V}as  an  LNs  x  1  vector  containing  all  signals  received 
by  all  L  sensors  when  the  emitter  is  in  grid  point  n  as 
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d„  =  r  x  .s' . 


(20) 


If  we  arrange  all  vectors  6n  for  n:l...N  as  the  columns  of  a  matrix  0  as 


0  =  [0,  02 


0N ] 


N  -ILN-xN  ? 


(21) 


and  then,  we  have 


with 


r  =  0  x  z  +  w 

II 

T 

r2  • 

r  T 

rL  J  LNsx  1 

Z  =  [z{ 

z2 

~\T 

•  ZN  J  Nx\  •> 

(22) 


where  r  is  the  vector  of  all  L  received  signals,  z  is  the  sparse  vector  of  z-values  assigned  to  each 
grid  point  and  w  is  the  noise.  Now,  we  can  solve  our  problem  by  fonning  a  BPIC  ( Basis  Pursuit 
with  Inequality  Constraints)  problem  [33]  as  following: 


(z  =  arginin^l 

I  s.t  ||0xz-r||2  <  s 

or  regularized  BPDN  ( Basis  Pursuit  Denoising)  problem  [33]  as: 

Z  =  arg  min  \\0  x  z  -  r ||2  +  A  1^1 


(23) 


(24) 


We  examined  the  performance  of  the  proposed  method  and  compared  the  results  using 
Monte-Carlo  computer  simulations  for  different  scenarios.  In  the  first  simulation,  we  assumed 
that  3  moving  sensors  receive  the  signal  from  one  stationary  emitter  placed  in  a  configuration  as 
shown  in  Figure  12  (the  location  of  the  emitter  has  been  chosen  randomly).  In  this  simulation, 
the  sampling  frequency  is  80  kHz  and  the  number  of  samples  is  equal  to  1024.  In  Figure  13,  we 
can  see  the  RMS  Error  vs.  SNR  (with  500  runs  for  each  SNR)  for  estimating  the  location  of  the 
emitter  in  ( x  -  y)  plane.  As  we  see,  the  proposed  method  has  better  perfonnance  compared  to 
DPD  and  Classic  methods. 
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Figure  12:  Placement  of  the  sensors  and  the  emitter  position  used  for  simulation. 


Figure  13:  Location  error  versus  SNR  for  single-emitter  single-path  case. 


One  of  the  challenging  topics  in  source  localization  problems  is  emitter  location  estimation  in 
the  presence  of  multipath  reflections.  We  evaluated  the  capability  of  the  proposed  method  in 
dealing  with  multipath  scenarios  using  Monte-Carlo  simulation.  In  this  simulation,  we  assumed 
that  4  moving  sensors  receive  the  signal  from  one  stationary  emitter  placed  in  a  configuration  as 
shown  in  Figure  14  (the  location  of  the  emitter  has  been  chosen  randomly).  Figure  15  shows  the 
RMS  Error  vs.  SNR  (with  500  runs  for  each  SNR)  for  estimating  the  location  of  the  emitter  in 
multipath  case.  The  following  plots  show  better  accuracy  of  the  proposed  method  over  DPD  and 
Classic  methods.  As  we  see,  none  of  the  DPD  and  Classic  methods  provide  unbiased  estimates 
for  higher  SNRs. 
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Figure  14:  Placement  of  the  sensors  and  the  emitter  position  used  for  simulation. 


Figure  15:  Location  error  versus  SNR  for  multipath  scenario. 


In  another  simulation,  we  evaluated  the  performance  of  the  proposed  method  and  compared 
the  results  to  DPD  and  Classic  method  for  multi-emitter  scenarios  when  we  have  more  than  one 
transmitter  in  the  range  of  interest  and  we  aim  to  estimate  the  location  of  all  emitters.  In  this 
simulation,  we  assumed  that  6  moving  sensors  receive  the  signal  from  two  stationary  emitters 
placed  in  the  configuration  as  shown  in  Figure  16  (the  location  of  the  emitters  has  been  chosen 
randomly).  Figure  17  shows  the  RMS  Error  vs.  SNR  (with  500  runs  for  each  SNR)  for  estimating 
the  location  of  two  emitters.  Comparing  the  RMS  Error  curves,  we  see  that  the  proposed  method 
achieves  much  better  accuracy  with  significantly  smaller  error  compared  to  DPD  and  Classic 
methods.  As  we  see,  similar  to  multipath  case,  both  DPD  and  Classic  methods  provide  biased 
estimates  at  high  SNRs. 
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Figure  16:  Placement  of  the  sensors  and  the  emitter  position  used  for  simulation. 


(a) 


(c) 


(b) 


(d) 


Figure  17:  Performance  for  locating  two  emitters,  (a),  (b)  1st  emitter;  (c),  (d)  2nd  emitter. 
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We  developed  a  one-stage  TDOA/FDOA  localization  method  based  on  spatial  sparsity  of 
emitters.  In  this  method,  we  assign  a  non-zero  number  to  each  one  of  the  grids  containing  an 
emitter  and  zero  to  the  rest  of  the  grid  points.  Thus,  the  vector  fonned  from  these  numbers  will 
be  a  sparse  unknown  vector  that  we  aim  to  estimate.  Since  each  element  of  this  vector 
corresponds  to  one  grid  point  in  (x,y)  plane,  we  can  estimate  the  location  of  emitters  by 
extracting  the  position  of  non-zero  elements  of  the  sparsest  vector  that  satisfy  the  TDOA/FDOA 
relationship  between  transmitted  signals  and  received  signals.  We  evaluated  the  performance  of 
the  proposed  method  using  Monte-Carlo  simulation.  Comparing  the  three  curves  in  each  plot  in 
Figure  13,  Figure  15  and  Figure  17  shows  that  the  proposed  method  has  better  perfonnance 
(especially  in  multi-path  and  multi-emitter  cases)  compared  to  direct  position  determination 
(DPD)  and  two-stage  Classic  localization  methods.  Simulation  results  show  that  contrary  to  DPD 
and  Classic  methods,  the  proposed  method  is  a  very  reliable  and  strong  tool  to  deal  with 
multipath  and  multi-emitter  scenarios. 


5.0  CONCLUSIONS 

For  several  decades  now  emitter  location  has  been  done  using  a  two-stage  approach  but 
recently  it  has  been  shown  that  a  single-stage  method  (called  DPD)  has  some  performance 
advantages.  However,  this  comes  at  the  cost  of  some  difficulties  that  hinder  easy 
implementation.  In  particular,  the  original  fonnulation  of  DPD  calls  for  the  transmission  of  all 
signals  to  a  single  sensor.  This  causes  several  problems:  (i)  communication  limitations  restrict 
transmission  of  all  data  to  a  single  sensor,  (ii)  all  computation  is  perfonned  at  a  single  sensor  - 
resulting  in  an  excessive  burden  for  one  sensor  and  (iii)  if  that  single  computing  sensor  is 
prevented  from  completing  the  location  (due  to  failure  or  enemy  action)  then  no  location  can  be 
made  available.  The  first  set  of  results  provided  here  address  these  issues  and  make  strides 
toward  a  feasible  implementation  of  the  DPD  method  -  thus  helping  to  make  DPD’s  improved 
perfonnance  more  readily  available  in  practice. 

However,  in  scenarios  of  multipath  and  multiple  co-channel  emitters  we  demonstrated  that 
even  though  the  DPD  method  does  improve  upon  the  classical  methods  it  still  leaves  a  need  for 
further  improvement.  We  showed  that  it  is  possible  to  implement  a  single-stage  method  that  is 
similar  to  DPD  but  that  also  exploits  the  sparsity  of  emitters  in  space.  The  simulation  results 
showed  that  this  method  significantly  improves  the  perfonnance  buy  at  least  2  orders  of 
magnitude!  This  is  a  tremendous  improvement  and  -  once  developed  to  be  feasible  in  the  real 
world  -  will  provide  accurate  locations  in  environments  previously  found  to  be  unworkable. 
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APPENDIX 


As  mentioned  above,  we  can  consider  the  CAF  as  an  image  and  we  can  apply  image 
compression  methods  to  compress  it.  However,  not  all  of  the  CAF  points  have  the  same 
importance  for  location  estimation,  thus  it  is  possible  to  assign  different  weights  to  different  CAF 
points  and  therefore  allocate  larger  number  of  data  bits  to  transmit  the  more  significant  area, 
which  contains  the  mainlobe  area. 

Price  and  Hofstetter  [12]  have  done  detailed  research  on  the  bounds  of  the  ambiguity 
function  volume  distribution.  Wilcox  [13]  showed  that  the  contour  of  ambiguity  function 
magnitude  close  to  the  peak  is  always  an  ellipse.  This  contour  can  be  formed  by  the  intersection 
of  the  mainlobe  magnitude  and  a  level  plane.  It  is  possible  to  find  the  approximate  equation  of 
this  ellipse  in  terms  of  signal  bandwidth,  signal  duration,  signal  energy  and  a  specific  level.  The 
width  of  this  ellipse  along  the  TDOA  axis  is  proportional  to  the  reciprocal  of  the  signal’s  rms 
bandwidth;  likewise,  the  width  of  the  ellipse  along  the  FDOA  axis  is  proportional  to  the 
reciprocal  of  signal’s  nns  duration  [14]: 


At  oc  (1  / Brms) 

Am  ex  (1  /Trms) 

where  Brms  is  rms  value  of  signal  bandwidth  and  Trms  is  rms  value  of  signal  duration. 

Thus,  it  is  possible  to  determine  the  approximate  significant  area  which  is  more  important  for 
the  purpose  of  location  estimation.  Note  that  it  is  always  possible  to  rotate  the  ambiguity  function 
when  it  is  tilted.  An  interesting  property  of  ambiguity  function  is  that  the  new  function  under  the 
transformation  that  rotates  the  r  —  m  plane  through  some  angle  9  will  be  another  ambiguity 
function.  This  new  ambiguity  function  is  corresponding  to  the  new  signals  which  are  related  to 
the  old  signals  and  the  angle  9  [14], [15]. 

Under  the  so-called  narrowband  approximation  the  lowpass  equivalent  (LPE)  model  of  the 
received  signal  will  be: 


Sr(t)  =  e  J'wcTd  e  J'wdts(t  —  Td )  (Al) 

where  s(t)  is  the  LPE  of  the  transmitted  signal,  wd  is  the  Doppler  and  rd  is  the  delay  for  the 
received  signal  [11].  Now,  suppose  that  two  sensors  Rxl  and  Rx2  receive  the  LPE  signals 
srl(t)  andsr2(t),  respectively.  The  ML  estimate  for  TDOA  and  FDOA  can  be  obtained  using 
the  magnitude  of  the  CAF: 

Ai2(r,m)=  C  srl(t)sr2*(t  -  r)ejcotdt  (A2) 

which  measures  the  correlation  between  srl  (t)  and  a  Doppler-shifted  by  co  and  delayed  by  r 
version  of  sr2(t). 

The  auto  ambiguity  function  (AAF)  (or  just  ambiguity  function  as  in  some  papers)  is  defined 
as: 
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AuuQr,  <w)  =  J_+”  u(t)u*(t  -  T )eja>tdt 


(A3) 


where  u(t)  can  be  the  LPE  signal.  In  fact,  auto  ambiguity  function  shows  the  correlation 
between  a  signal  and  a  Doppler-shifted  by  u>  and  delayed  by  r  version  of  itself.  It  is  straight 
forward  to  show  that  AAF  has  a  kind  of  symmetry  around  the  origin  [13],  [14],  [16]. 

Auu(.-t,-<o)  =  A*uu(  (A4) 

|Auu(  A  tu)  |  \Auu(t,0))\ 

where  A^^t,  u>)  is  the  complex  conjugate  of  AAF  and  \Auu(t,  u>)  |  is  the  magnitude  of  AAF. 

It  is  also  simple  to  prove  a  similar  property  for  CAF  [14], [16]: 

Auv(-t,-o))  =  A*vu( t,u>)  ejT0J 

A  tu)  |  |y4vu(T,  U))| 

where  Auv(r,  u>)  is  the  CAF  between  signal  u(t)  and  v(t).  However,  Avu(t,co)  is  the  CAF 
between  arbitrary  signals  v(t )  and  it(t),  not  specifically  related  by  delay  and  Doppler  to  a  single 
transmitted  signal.  To  develop  a  result  that  we  can  exploit  for  our  purpose  we  explore  a  similar 
result  for  the  case  when  the  signals  are  received  from  a  transmitter.  Then  (Al)  gives 

it(t)  =  e~J(x>cTl  e~Jlx>lt  s(t  —  rx) 

v(t)  —  e~J0)cT2  e~J0)2t  s(t  —  r2 ) 

where  and  r2  are  the  time  delays  and  oq  and  o)2  are  the  Doppler  shifts  for  the  first  and  second 
received  signals.  Now,  we  can  write  one  of  them  in  terms  of  the  other  one, 

v(t)  —  it(t  +  r p)e-,WcTp  eJ(x>lTp  (A5) 

where  rp  =  (rx  —  r2)  is  the  TDOA  and  cop  —  (uq  —  u>2)  is  the  FDOA. 

Auv(j,  <u)  =  /_+“  u(t)v*(t  -  r )eja)tdt  S 

=  ej*>P*-J*>cTp-J<oiTv  J_+“w(t)u*(t  -  (-Tp  +  T))  e^-^dt 

=  [etwpT-j«crp-ja)iTpj  -  Tp>  U)  -  (Dp)  (A6) 

Thus,  the  CAF  is  rewritten  in  terms  of  AAF.  Then,  by  replacing  the  r  by  (t  +  rp)  and  <u 
by  (<u  —  <Up),  the  following  equations  are  concluded, 
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(A7) 


Auv( r  +  r p,o)  +  cop)  =  eJ'wp(T+Tp)  ja>cTp  ja>lTp  Auu(t,o) ) 

<d)  =  e;'wp(T+Tp)_;'WcTP_;'"lTP  i4yV(r  +  rp,  <d  +  op)  (A8) 

Now,  by  negating  the  r  and  <d  in  (A7),  we  have: 


Auv  (  f  T  ft)  T  ft)p  ) 

=  [e-/wp(-T+Tp)-J&)cTp-;&)iTpj  Auu(  —  T,  —(D) 


(^4) 


—  JeJ'<Wp(-T+Tp)-/ 


7COcTp  JOJ^Tp  p 


e  J0)T]A*UU(  T,a>) 


(A  8) 


—  T+Tp)  i(^cTp  j^^Tp  Q  j(x) T  gj(Op{r+T-^)  jk)cTp  JM^Tp  -|-  -|-  (x)p^ 


=  [eU2*>P*p-j2<*CTp-j2<OlTp)  e-jC0 T]  ^(T  +  ^  +  <yp) 

=  e-t(WT+^)^v(T  +  Tp)W  +  <up) 

where  /?  is  defined  as  (2copTp  —  2cocrp  —  2u)1rp')  and  finally, 

Auv(-t  +  tp>  -ft)  +  ft)p)  =  e~j(ojr+f3)  A*uv(t  +  tp,u)  +  op) 
and 

I  Auv(  T  +  Tp,  (D  +  <Dp)  I  I  Auv(t  +  Tp,  (D  +  <Dp)  | 


(A9) 

(AlO) 


which  is  the  symmetry  property  we  can  exploit  for  data  compression.  This  result  provides  a  kind 
of  symmetry  of  the  CAF  around  the  point  (rp,  cop)  or  the  peak  of  CAF  magnitude. 

Now,  it  is  possible  to  exploit  this  property  in  data  compression.  In  practice,  the  received 
signals  u(t )  and  v(t)  are  the  delayed  and  Doppler-shifted  version  of  transmitted  signal  plus 
noise.  This  noise  perturbs  the  CAF  a  little  bit  from  the  perfect  symmetry. 

Thus,  we  rewrite  (AlO)  as, 

I  AUi;(  T  +  Tp,  (D  +  <Dp)  |  |  Auv(t  +  Tp,  CO  +  <Dp)  |  +  E  (A1 1) 

where  E  can  be  the  error  from  perfect  symmetry  which  is  a  negligible  value.  Thus,  using  the 
symmetry  property,  it  is  possible  to  extract  the  entire  CAF  magnitude  by  transmission  of  only 
half  of  the  CAF  magnitude  plus  the  small  residual  amount  of  E.  In  this  scheme  we  apply  the 
EZW  data  compression  method  on  only  half  of  CAF  as  well  as  on  E. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS  AND  ACRONYMS 


AAF 

auto  ambiguity  function 

ADC 

analog-to-digital  converter 

AESA 

active  electronically  scanned  array 

BPDN 

basis  pursuit  denoising 

BPIC 

basis  pursuit  with  inequality  constraints 

BPSK 

binary  phase  shift  keying 

CAF 

cross-ambiguity  function 

CAM 

cross  ambiguity  matrix 

COMINT 

communication  intelligence 

COTS 

commercial-off-the-shelf 

CRLB 

Cramer-Rao  lower  bound 

DAWG 

digital  arbitrary  waveform  generator 

DPD 

direct  position  detennination 

ELINT 

electronics  intelligence 

EZW 

embedded  zerotree  wavelet 

FDOA 

frequency-difference-of-arrival 

GLR 

generalized  likelihood  ratio 

LPE 

lowpass  equivalent 

LPI 

low-probability  of  intercept 

ML 

maximum  likelihood 

RF 

radio  frequency 

SNR 

signal-to-noise  ratio 

SVD 

singular  value  decomposition 

TDOA 

time-difference-of-arrival 
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